Data Import

data <- readRDS("very_low_birthweight.RDS")

Task 1

Загрузите датасет very_low_birthweight.RDS (лежит в папке домашнего задания). Это данные о 671 младенце с очень низкой массой тела (<1600 грамм), собранные в Duke University Medical Center доктором Майклом О’Ши c 1981 по 1987 г.
Описание переменных см. здесь. Переменными исхода являются колонки ‘dead’, а также время от рождения до смерти или выписки (выводятся из ‘birth’ и ‘exit’ 7 пациентов были выписаны до рождения). Сделайте копию датасета, в которой удалите колонки с количеством пропусков больше 100, а затем удалите все строки с пропусками.

data %>% 
  head()
##    birth   exit hospstay    lowph pltct  race  bwt gest        inout twn lol
## 1 81.511 81.604       34       NA   100 white 1250   35 born at Duke   0  NA
## 2 81.514 81.539        9 7.250000   244 white 1370   32 born at Duke   0  NA
## 3 81.552 81.552       -2 7.059998   114 black  620   23 born at Duke   0  NA
## 4 81.558 81.667       40 7.250000   182 black 1480   32 born at Duke   0  NA
## 5 81.593 81.599        2 6.969997    54 black  925   28 born at Duke   0  NA
## 6 81.602 81.771       62 7.189999    NA white  940   28 born at Duke   0  NA
##   magsulf meth toc  delivery apg1 vent pneumo pda cld      pvh      ivh    ipe
## 1      NA    0   0 abdominal    8    0      0   0   0     <NA>     <NA>   <NA>
## 2      NA    1   0 abdominal    7    0      0   0   0     <NA>     <NA>   <NA>
## 3      NA    0   1   vaginal    1    1      0   0  NA     <NA>     <NA>   <NA>
## 4      NA    1   0   vaginal    8    0      0   0   0     <NA>     <NA>   <NA>
## 5      NA    0   0 abdominal    5    1      1   0   0 definite definite   <NA>
## 6      NA    1   0 abdominal    8    1      0   0   0   absent   absent absent
##       year    sex dead
## 1 81.51196 female    0
## 2 81.51471 female    0
## 3 81.55304 female    1
## 4 81.55847   male    0
## 5 81.59406 female    1
## 6 81.60229 female    0
cols_to_select <- data.frame(NA_count = colSums(is.na(data))) %>% filter(NA_count <= 100) %>% row.names()

data %>% 
  select(all_of(cols_to_select)) %>% 
  na.omit() %>% 
  write_csv("filtered_dataset.csv")

Task 2

Постройте графики плотности распределения для числовых переменных. Удалите выбросы, если таковые имеются. Преобразуйте категориальные переменные в факторы. Для любых двух числовых переменных раскрасьте график по переменной ‘inout’.

data <- data %>% 
  mutate(across(where(is.character), as.numeric), 
         across(where(is.factor), as.factor),
         across(where(~ is.numeric(.) & n_distinct(., na.rm = TRUE) == 2), as.factor))

data %>% 
  select_if(is.numeric) %>% 
  pivot_longer(cols = everything(), names_to = "variable", values_to = "value") %>% 
  ggplot(aes(x = value)) +
  geom_histogram(bins = 10, fill = "skyblue", color = "black") +
  facet_wrap(~variable, scales = "free_x") + 
  theme_minimal() +
  theme(axis.text.x = element_text(size=7))
## Warning: Removed 657 rows containing non-finite outside the scale range
## (`stat_bin()`).

dfSummary(data)
## Data Frame Summary  
## data  
## Dimensions: 671 x 26  
## Duplicates: 0  
## 
## --------------------------------------------------------------------------------------------------------------
## No   Variable    Stats / Values               Freqs (% of Valid)    Graph                 Valid      Missing  
## ---- ----------- ---------------------------- --------------------- --------------------- ---------- ---------
## 1    birth       Mean (sd) : 84.8 (1.6)       542 distinct values             : :         650        21       
##      [numeric]   min < med < max:                                         :   : : : . :   (96.9%)    (3.1%)   
##                  81.5 < 84.9 < 87.5                                   : : : : : : : : :                       
##                  IQR (CV) : 2.6 (0)                                 : : : : : : : : : :                       
##                                                                     : : : : : : : : : :                       
## 
## 2    exit        Mean (sd) : 84.8 (1.8)       539 distinct values             :           640        31       
##      [numeric]   min < med < max:                                             : .         (95.4%)    (4.6%)   
##                  68.5 < 85 < 96.9                                             : :                             
##                  IQR (CV) : 2.6 (0)                                           : :                             
##                                                                             : : :                             
## 
## 3    hospstay    Mean (sd) : 40.4 (304.8)     153 distinct values               :         640        31       
##      [integer]   min < med < max:                                               :         (95.4%)    (4.6%)   
##                  -6574 < 37 < 3668                                              :                             
##                  IQR (CV) : 46 (7.6)                                            :                             
##                                                                                 :                             
## 
## 4    lowph       Mean (sd) : 7.2 (0.1)        74 distinct values                : .       609        62       
##      [numeric]   min < med < max:                                               : :       (90.8%)    (9.2%)   
##                  6.5 < 7.2 < 7.5                                              . : :                           
##                  IQR (CV) : 0.2 (0)                                         . : : : :                         
##                                                                           . : : : : : .                       
## 
## 5    pltct       Mean (sd) : 201.6 (80.5)     266 distinct values         :               601        70       
##      [integer]   min < med < max:                                       . : .             (89.6%)    (10.4%)  
##                  16 < 202 < 571                                       . : : :                                 
##                  IQR (CV) : 109 (0.4)                                 : : : : .                               
##                                                                     : : : : : : .                             
## 
## 6    race        1. black                     369 (57.1%)           IIIIIIIIIII           646        25       
##      [factor]    2. native American            16 ( 2.5%)                                 (96.3%)    (3.7%)   
##                  3. oriental                    4 ( 0.6%)                                                     
##                  4. white                     257 (39.8%)           IIIIIII                                   
## 
## 7    bwt         Mean (sd) : 1093.9 (265.2)   142 distinct values               . :       669        2        
##      [integer]   min < med < max:                                           . : : : :     (99.7%)    (0.3%)   
##                  400 < 1120 < 1580                                        : : : : : :                         
##                  IQR (CV) : 410 (0.2)                                 . : : : : : : : :                       
##                                                                     . : : : : : : : : :                       
## 
## 8    gest        Mean (sd) : 28.9 (2.5)       24 distinct values        : :               667        4        
##      [numeric]   min < med < max:                                       : : :             (99.4%)    (0.6%)   
##                  22 < 29 < 40                                         . : : :                                 
##                  IQR (CV) : 4 (0.1)                                   : : : :                                 
##                                                                     . : : : : :                               
## 
## 9    inout       1. born at Duke              547 (81.9%)           IIIIIIIIIIIIIIII      668        3        
##      [factor]    2. transported               121 (18.1%)           III                   (99.6%)    (0.4%)   
## 
## 10   twn         1. 0                         516 (79.3%)           IIIIIIIIIIIIIII       651        20       
##      [factor]    2. 1                         135 (20.7%)           IIII                  (97.0%)    (3.0%)   
## 
## 11   lol         Mean (sd) : 8.4 (19.3)       40 distinct values    :                     290        381      
##      [integer]   min < med < max:                                   :                     (43.2%)    (56.8%)  
##                  0 < 3.5 < 192                                      :                                         
##                  IQR (CV) : 9 (2.3)                                 :                                         
##                                                                     : .                                       
## 
## 12   magsulf     1. 0                         367 (86.6%)           IIIIIIIIIIIIIIIII     424        247      
##      [factor]    2. 1                          57 (13.4%)           II                    (63.2%)    (36.8%)  
## 
## 13   meth        1. 0                         318 (56.3%)           IIIIIIIIIII           565        106      
##      [factor]    2. 1                         247 (43.7%)           IIIIIIII              (84.2%)    (15.8%)  
## 
## 14   toc         1. 0                         438 (77.5%)           IIIIIIIIIIIIIII       565        106      
##      [factor]    2. 1                         127 (22.5%)           IIII                  (84.2%)    (15.8%)  
## 
## 15   delivery    1. abdominal                 314 (48.4%)           IIIIIIIII             649        22       
##      [factor]    2. vaginal                   335 (51.6%)           IIIIIIIIII            (96.7%)    (3.3%)   
## 
## 16   apg1        Mean (sd) : 4.9 (2.6)        0 :   5 ( 0.8%)                             637        34       
##      [integer]   min < med < max:             1 :  92 (14.4%)       II                    (94.9%)    (5.1%)   
##                  0 < 5 < 9                    2 :  69 (10.8%)       II                                        
##                  IQR (CV) : 5 (0.5)           3 :  57 ( 8.9%)       I                                         
##                                               4 :  45 ( 7.1%)       I                                         
##                                               5 :  69 (10.8%)       II                                        
##                                               6 :  80 (12.6%)       II                                        
##                                               7 :  80 (12.6%)       II                                        
##                                               8 : 103 (16.2%)       III                                       
##                                               9 :  37 ( 5.8%)       I                                         
## 
## 17   vent        1. 0                         269 (42.0%)           IIIIIIII              641        30       
##      [factor]    2. 1                         372 (58.0%)           IIIIIIIIIII           (95.5%)    (4.5%)   
## 
## 18   pneumo      1. 0                         518 (80.3%)           IIIIIIIIIIIIIIII      645        26       
##      [factor]    2. 1                         127 (19.7%)           III                   (96.1%)    (3.9%)   
## 
## 19   pda         1. 0                         508 (79.1%)           IIIIIIIIIIIIIII       642        29       
##      [factor]    2. 1                         134 (20.9%)           IIII                  (95.7%)    (4.3%)   
## 
## 20   cld         1. 0                         442 (73.1%)           IIIIIIIIIIIIII        605        66       
##      [factor]    2. 1                         163 (26.9%)           IIIII                 (90.2%)    (9.8%)   
## 
## 21   pvh         1. absent                    360 (68.4%)           IIIIIIIIIIIII         526        145      
##      [factor]    2. definite                  125 (23.8%)           IIII                  (78.4%)    (21.6%)  
##                  3. possible                   41 ( 7.8%)           I                                         
## 
## 22   ivh         1. absent                    442 (83.9%)           IIIIIIIIIIIIIIII      527        144      
##      [factor]    2. definite                   75 (14.2%)           II                    (78.5%)    (21.5%)  
##                  3. possible                   10 ( 1.9%)                                                     
## 
## 23   ipe         1. absent                    472 (89.6%)           IIIIIIIIIIIIIIIII     527        144      
##      [factor]    2. definite                   38 ( 7.2%)           I                     (78.5%)    (21.5%)  
##                  3. possible                   17 ( 3.2%)                                                     
## 
## 24   year        Mean (sd) : 84.8 (1.6)       542 distinct values             : :         650        21       
##      [numeric]   min < med < max:                                         :   : : : . :   (96.9%)    (3.1%)   
##                  81.5 < 84.9 < 87.5                                   : : : : : : : : :                       
##                  IQR (CV) : 2.6 (0)                                 : : : : : : : : : :                       
##                                                                     : : : : : : : : : :                       
## 
## 25   sex         1. female                    320 (49.2%)           IIIIIIIII             650        21       
##      [factor]    2. male                      330 (50.8%)           IIIIIIIIII            (96.9%)    (3.1%)   
## 
## 26   dead        1. 0                         527 (78.5%)           IIIIIIIIIIIIIII       671        0        
##      [factor]    2. 1                         144 (21.5%)           IIII                  (100.0%)   (0.0%)   
## --------------------------------------------------------------------------------------------------------------

Что-то странное происходит в hospstay, предполагаю, что младенцы не умеют путешествовать во времени. В этих (exit, pltct, lol) переменных так же есть экстремальные значения, природа которых мне неизвестна, поэтому удалю их по правилу трех IQR.

NB: после удаления аутлаеров по hospstay исчезли аутлаеры у exit

data <- data %>% 
  mutate(ID = c(1:nrow(data))) %>% 
  select(ID, everything()) 
data <- data %>% 
  filter(hospstay > 0 | is.na(hospstay), 
         hospstay < 3 * IQR(hospstay, na.rm = T) | is.na(hospstay),
         pltct < 3 * IQR(pltct, na.rm = T) | is.na(pltct),
         lol < 3 * IQR(lol, na.rm = T) | is.na(lol))
data %>% 
  dfSummary()
## Data Frame Summary  
## Dimensions: 569 x 27  
## Duplicates: 0  
## 
## ---------------------------------------------------------------------------------------------------------------
## No   Variable    Stats / Values               Freqs (% of Valid)    Graph                  Valid      Missing  
## ---- ----------- ---------------------------- --------------------- ---------------------- ---------- ---------
## 1    ID          Mean (sd) : 332.9 (193.8)    569 distinct values   . : : . : . . . . .    569        0        
##      [integer]   min < med < max:                                   : : : : : : : : : :    (100.0%)   (0.0%)   
##                  1 < 330 < 671                                      : : : : : : : : : :                        
##                  IQR (CV) : 336 (0.6)                               : : : : : : : : : :                        
##                                                                     : : : : : : : : : :                        
## 
## 2    birth       Mean (sd) : 84.7 (1.6)       460 distinct values             : .          548        21       
##      [numeric]   min < med < max:                                         :   : : :   .    (96.3%)    (3.7%)   
##                  81.5 < 84.9 < 87.5                                   : : : : : : : : :                        
##                  IQR (CV) : 2.5 (0)                                 : : : : : : : : : :                        
##                                                                     : : : : : : : : : :                        
## 
## 3    exit        Mean (sd) : 84.8 (1.6)       460 distinct values             :            538        31       
##      [numeric]   min < med < max:                                         : . : : :        (94.6%)    (5.4%)   
##                  81.5 < 84.9 < 87.8                                   : : : : : : : : .                        
##                  IQR (CV) : 2.6 (0)                                 . : : : : : : : : :                        
##                                                                     : : : : : : : : : :                        
## 
## 4    hospstay    Mean (sd) : 41.6 (30.6)      117 distinct values   :                      538        31       
##      [integer]   min < med < max:                                   : : : .                (94.6%)    (5.4%)   
##                  1 < 36 < 137                                       : : : : .                                  
##                  IQR (CV) : 43 (0.7)                                : : : : : .                                
##                                                                     : : : : : : : . . .                        
## 
## 5    lowph       Mean (sd) : 7.2 (0.1)        71 distinct values                : .        516        53       
##      [numeric]   min < med < max:                                               : :        (90.7%)    (9.3%)   
##                  6.5 < 7.2 < 7.5                                              . : :                            
##                  IQR (CV) : 0.2 (0)                                         . : : : :                          
##                                                                           . : : : : : .                        
## 
## 6    pltct       Mean (sd) : 192.1 (68.1)     222 distinct values           :              509        60       
##      [integer]   min < med < max:                                         : :              (89.5%)    (10.5%)  
##                  16 < 199 < 324                                         : : : .                                
##                  IQR (CV) : 99 (0.4)                                  . : : : :                                
##                                                                     . : : : : : :                              
## 
## 7    race        1. black                     303 (55.7%)           IIIIIIIIIII            544        25       
##      [factor]    2. native American            16 ( 2.9%)                                  (95.6%)    (4.4%)   
##                  3. oriental                    3 ( 0.6%)                                                      
##                  4. white                     222 (40.8%)           IIIIIIII                                   
## 
## 8    bwt         Mean (sd) : 1110.1 (258.9)   134 distinct values                 :        567        2        
##      [integer]   min < med < max:                                             . : : :      (99.6%)    (0.4%)   
##                  400 < 1150 < 1580                                        . : : : : :                          
##                  IQR (CV) : 400 (0.2)                                   . : : : : : : .                        
##                                                                     . : : : : : : : : :                        
## 
## 9    gest        Mean (sd) : 29.1 (2.5)       23 distinct values        . :                566        3        
##      [numeric]   min < med < max:                                       : : :              (99.5%)    (0.5%)   
##                  23 < 29 < 40                                           : : :                                  
##                  IQR (CV) : 4 (0.1)                                   : : : :                                  
##                                                                     . : : : : :                                
## 
## 10   inout       1. born at Duke              462 (81.5%)           IIIIIIIIIIIIIIII       567        2        
##      [factor]    2. transported               105 (18.5%)           III                    (99.6%)    (0.4%)   
## 
## 11   twn         1. 0                         432 (78.5%)           IIIIIIIIIIIIIII        550        19       
##      [factor]    2. 1                         118 (21.5%)           IIII                   (96.7%)    (3.3%)   
## 
## 12   lol         Mean (sd) : 5.2 (6.3)        26 distinct values    :                      242        327      
##      [integer]   min < med < max:                                   :                      (42.5%)    (57.5%)  
##                  0 < 3 < 26                                         :                                          
##                  IQR (CV) : 8 (1.2)                                 : .                                        
##                                                                     : : . . .                                  
## 
## 13   magsulf     1. 0                         306 (85.2%)           IIIIIIIIIIIIIIIII      359        210      
##      [factor]    2. 1                          53 (14.8%)           II                     (63.1%)    (36.9%)  
## 
## 14   meth        1. 0                         263 (55.4%)           IIIIIIIIIII            475        94       
##      [factor]    2. 1                         212 (44.6%)           IIIIIIII               (83.5%)    (16.5%)  
## 
## 15   toc         1. 0                         372 (78.5%)           IIIIIIIIIIIIIII        474        95       
##      [factor]    2. 1                         102 (21.5%)           IIII                   (83.3%)    (16.7%)  
## 
## 16   delivery    1. abdominal                 274 (50.1%)           IIIIIIIIII             547        22       
##      [factor]    2. vaginal                   273 (49.9%)           IIIIIIIII              (96.1%)    (3.9%)   
## 
## 17   apg1        Mean (sd) : 5 (2.6)          0 :  3 ( 0.6%)                               536        33       
##      [integer]   min < med < max:             1 : 74 (13.8%)        II                     (94.2%)    (5.8%)   
##                  0 < 5 < 9                    2 : 53 ( 9.9%)        I                                          
##                  IQR (CV) : 4 (0.5)           3 : 46 ( 8.6%)        I                                          
##                                               4 : 35 ( 6.5%)        I                                          
##                                               5 : 61 (11.4%)        II                                         
##                                               6 : 68 (12.7%)        II                                         
##                                               7 : 71 (13.2%)        II                                         
##                                               8 : 96 (17.9%)        III                                        
##                                               9 : 29 ( 5.4%)        I                                          
## 
## 18   vent        1. 0                         236 (43.5%)           IIIIIIII               543        26       
##      [factor]    2. 1                         307 (56.5%)           IIIIIIIIIII            (95.4%)    (4.6%)   
## 
## 19   pneumo      1. 0                         441 (81.1%)           IIIIIIIIIIIIIIII       544        25       
##      [factor]    2. 1                         103 (18.9%)           III                    (95.6%)    (4.4%)   
## 
## 20   pda         1. 0                         441 (81.1%)           IIIIIIIIIIIIIIII       544        25       
##      [factor]    2. 1                         103 (18.9%)           III                    (95.6%)    (4.4%)   
## 
## 21   cld         1. 0                         393 (76.2%)           IIIIIIIIIIIIIII        516        53       
##      [factor]    2. 1                         123 (23.8%)           IIII                   (90.7%)    (9.3%)   
## 
## 22   pvh         1. absent                    307 (69.3%)           IIIIIIIIIIIII          443        126      
##      [factor]    2. definite                  100 (22.6%)           IIII                   (77.9%)    (22.1%)  
##                  3. possible                   36 ( 8.1%)           I                                          
## 
## 23   ivh         1. absent                    378 (85.3%)           IIIIIIIIIIIIIIIII      443        126      
##      [factor]    2. definite                   57 (12.9%)           II                     (77.9%)    (22.1%)  
##                  3. possible                    8 ( 1.8%)                                                      
## 
## 24   ipe         1. absent                    399 (90.1%)           IIIIIIIIIIIIIIIIII     443        126      
##      [factor]    2. definite                   30 ( 6.8%)           I                      (77.9%)    (22.1%)  
##                  3. possible                   14 ( 3.2%)                                                      
## 
## 25   year        Mean (sd) : 84.7 (1.6)       460 distinct values             : .          548        21       
##      [numeric]   min < med < max:                                         :   : : :   .    (96.3%)    (3.7%)   
##                  81.5 < 84.9 < 87.5                                   : . : : : : : : :                        
##                  IQR (CV) : 2.5 (0)                                 . : : : : : : : : :                        
##                                                                     : : : : : : : : : :                        
## 
## 26   sex         1. female                    277 (50.5%)           IIIIIIIIII             548        21       
##      [factor]    2. male                      271 (49.5%)           IIIIIIIII              (96.3%)    (3.7%)   
## 
## 27   dead        1. 0                         462 (81.2%)           IIIIIIIIIIIIIIII       569        0        
##      [factor]    2. 1                         107 (18.8%)           III                    (100.0%)   (0.0%)   
## ---------------------------------------------------------------------------------------------------------------
data %>% 
  select_if(is.numeric) %>% 
  select(-ID) %>% 
  pivot_longer(cols = everything(), names_to = "variable", values_to = "value") %>% 
  ggplot(aes(x = value)) +
  geom_histogram(bins = 10, fill = "skyblue", color = "black") +
  facet_wrap(~variable, scales = "free_x") + 
  theme_minimal() +
  theme(axis.text.x = element_text(size=7))
## Warning: Removed 582 rows containing non-finite outside the scale range
## (`stat_bin()`).

Теперь выглядит приятнее :)

Task 3

Проведите тест на сравнение значений колонки ‘lowph’ между группами в переменной inout. Вид статистического теста определите самостоятельно. Визуализируйте результат через библиотеку ‘rstatix’. Как бы вы интерпретировали результат, если бы знали, что более низкое значение lowph ассоциировано с более низкой выживаемостью?

data %>% group_by(inout) %>% select(inout, lowph) %>% dfSummary()
## Data Frame Summary  
## Group: inout = born at Duke  
## Dimensions: 462 x 2  
## Duplicates: 396  
## 
## -------------------------------------------------------------------------------------------------------
## No   Variable    Stats / Values          Freqs (% of Valid)   Graph                 Valid     Missing  
## ---- ----------- ----------------------- -------------------- --------------------- --------- ---------
## 2    lowph       Mean (sd) : 7.2 (0.1)   65 distinct values               : :       420       42       
##      [numeric]   min < med < max:                                         : :       (90.9%)   (9.1%)   
##                  6.5 < 7.2 < 7.5                                        . : : .                        
##                  IQR (CV) : 0.2 (0)                                     : : : :                        
##                                                                     . : : : : : .                      
## -------------------------------------------------------------------------------------------------------
## 
## Group: inout = transported  
## Dimensions: 105 x 2  
## Duplicates: 60  
## 
## ---------------------------------------------------------------------------------------------------
## No   Variable    Stats / Values          Freqs (% of Valid)   Graph             Valid     Missing  
## ---- ----------- ----------------------- -------------------- ----------------- --------- ---------
## 2    lowph       Mean (sd) : 7.1 (0.2)   44 distinct values             :       96        9        
##      [numeric]   min < med < max:                                       : .     (91.4%)   (8.6%)   
##                  6.7 < 7.2 < 7.5                                        : :                        
##                  IQR (CV) : 0.2 (0)                                 : : : : .                      
##                                                               . . . : : : : :                      
## ---------------------------------------------------------------------------------------------------
## 
## Group: inout = NA  
## Dimensions: 2 x 2  
## Duplicates: 1  
## 
## ----------------------------------------------------------------------------------
## No   Variable    Stats / Values   Freqs (% of Valid)   Graph   Valid    Missing   
## ---- ----------- ---------------- -------------------- ------- -------- ----------
## 2    lowph       All NA's                                      0        2         
##      [numeric]                                                 (0.0%)   (100.0%)  
## ----------------------------------------------------------------------------------
data %>% 
  t_test(lowph ~ inout, detailed = T)
## Warning: The statistic is based on a difference or ratio; by default, for
## difference-based statistics, the explanatory variable is subtracted in the
## order "born at Duke" - "transported", or divided in the order "born at Duke" /
## "transported" for ratio-based statistics. To specify this order yourself,
## supply `order = c("born at Duke", "transported")`.
## # A tibble: 1 × 7
##   statistic  t_df     p_value alternative estimate lower_ci upper_ci
##       <dbl> <dbl>       <dbl> <chr>          <dbl>    <dbl>    <dbl>
## 1      5.20  130. 0.000000744 two.sided     0.0881   0.0546    0.122

У транспортированных пациентов значение рН меньше, чем у рожденных в Дюке (MD=0.088, CI 95%: [0.055; 0.122], p < \(10^{-6}\)). Поскольку мы знаем, что более низкое значение рН ассоциировано с более высокой смертностью, можно предположить, что транспортировка пациентов может быть также связана с более высокой смертностью.

Task 4

Сделайте новый датафрейм, в котором оставьте только континуальные или ранговые данные, кроме ‘birth’, ‘year’ и ‘exit’. Сделайте корреляционный анализ этих данных. Постройте два любых типа графиков для визуализации корреляций.

corr_data <- data %>% 
  select(where(is.numeric), -c(ID, birth, year, exit)) %>% 
  psych::corr.test(method = "spearman", adjust="BH")

cor_melted <- reshape2::melt(corr_data$r)
p_melted <- reshape2::melt(corr_data$p)

cor_plot_data <- merge(cor_melted, p_melted, by = c("Var1", "Var2"))
colnames(cor_plot_data) <- c("Var1", "Var2", "Correlation", "Adjusted_P")

ggplot(cor_plot_data, aes(Var1, Var2)) +
  geom_tile(aes(fill = Correlation), color = "white") +
  scale_fill_gradient2(low = "#433CE1", high = "#E13C59", mid = "grey90", 
                       midpoint = 0, limit = c(-1, 1), name = "Correlation") +
  geom_text(aes(label = ifelse(Adjusted_P < 0.05, sprintf("%.2f", Adjusted_P), "")), 
             color = "black", size = 4) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 12), 
        axis.text.y = element_text(size = 12)) +
  labs(x = "", y = "")

data %>% 
  select(where(is.numeric), -c(ID, birth, year, exit)) %>%
  ggpairs(progress = FALSE) +
  theme_bw() +
  theme(axis.text = element_text(size=6))

Task 5

Постройте иерархическую кластеризацию на этом датафрейме.

rownames(data) <- data$ID

res <- data %>% 
  select(where(is.numeric), -c(ID, birth, year, exit)) %>%
  na.omit() %>% 
  scale() %>% 
  dist(method = "euclidean") %>% 
  hclust(method = "ward.D2") 
res %>% 
  fviz_dend(cex = 0.2, 
            k = 3, 
            k_colors = "jco")
## Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
## of ggplot2 3.3.4.
## ℹ The deprecated feature was likely used in the factoextra package.
##   Please report the issue at <https://github.com/kassambara/factoextra/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Task 6

Сделайте одновременный график heatmap и иерархической кластеризации. Интерпретируйте результат.

data %>% 
  select(where(is.numeric), -c(ID, birth, year, exit)) %>%
  na.omit() %>% 
  scale() %>% 
  pheatmap::pheatmap(color=colorRampPalette(c("#433CE1", "grey90", "#E13C59"))(50), 
                     fontsize_row = 0.00001, fontsize_col = 7, angle_col = 0)

data %>% 
  select(where(is.numeric), -c(ID, birth, year, exit)) %>%
  na.omit() %>% 
  scale() %>% 
  pheatmap::pheatmap(color=colorRampPalette(c("#433CE1", "grey90", "#E13C59"))(50), 
                     fontsize_row = 0.00001, fontsize_col = 7, angle_col = 0, kmeans_k=4)

Визульно, можно выделить на хитмапе 1 4 кластера, чтобы их проинтерпретировать, можно усреднить значения в кластере, как на хитмапе 2. Первые 2 кластера довольно сильно похожи друг на друга, за исключением того, что пациенты из 1 кластера имели более выскойи уровень тромбоцитов и скор APG в отличие от пациентов из кластера 1. Пациенты из кластера 3 в отличие от остальных имеют наиболее длительный срок нахождения в больнице. Отличительными чертами пациентов из кластера 4 является самый низкий уровень рН, масса тела при рождении, а также минимальные сроки беременности.

Task 7 & 8

Проведите PCA анализ на этих данных. Проинтерпретируйте результат. Нужно ли применять шкалирование для этих данных перед проведением PCA?

Постройте biplot график для PCA. Раскрасьте его по значению колонки ‘dead’.

pca.data <- data %>% 
  select(where(is.numeric), -c(ID, birth, year, exit)) %>%
  na.omit() %>% 
  scale() %>% 
  prcomp()

fviz_eig(pca.data, 
         addlabels = T, 
         ylim = c(0, 50))

В нашем случае необходимо проводить шкалирование данных, так как все переменные имеют разные единицы измерения и разный размах. Первые две компоненты довольно хорошо объясняют данные, объясняя около 55% вариабельности.

pca.data$x %>% 
  ggplot() +
  geom_point(aes(x = PC1, y = PC2)) +
  theme_bw()

Вклад каждой переменной в компоненту:

pcs <- reshape2::melt(pca.data$rotation)


ggplot(pcs, aes(Var1, Var2)) +
  geom_tile(aes(fill = value), color = "white") +
  scale_fill_gradient2(low = "#433CE1", high = "#E13C59", mid = "grey90", 
                       midpoint = 0, limit = c(-1, 1)) +
  geom_text(aes(label = sprintf("%.2f", value)), 
             color = "black", size = 4) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 12), 
        axis.text.y = element_text(size = 12)) +
  labs(x = "", y = "")

В первую компоненту основной положительный вклад имеет время, проведенное в госпитале, и негативная масса тела при рождении и время беременности (что похоже на последний кластер в анализе выше).

Во вторую компоненту основной вклад вносит продолжительность родов. В третий уровень тромбоцитов и APG скор в большей степени и в меньшей степени продолжительность нахождения в больнице.

deaths <- data %>% filter(ID %in% rownames(pca.data$x)) %>% select(dead)

levels(deaths$dead) <- c("Alive", "Dead")

ggbiplot(pca.data, choices=1:2,
         scale=1, alpha = 0.5, groups = as.factor(deaths$dead)) + 
  theme_bw()

Task 9

Переведите последний график в ‘plotly’. При наведении на точку нужно, чтобы отображалось id пациента.

explained_variance <- pca.data$sdev
explained_variance_pct <- round(100 * explained_variance ** 2 / sum(explained_variance ** 2), 2)
comp <- pca.data$rotation
loadings <- comp
for (i in seq(explained_variance)){
  loadings[,i] <- comp[,i] * explained_variance[i]
}

features <- rownames(pca.data$rotation)

fig <- plot_ly(
  data = as.data.frame(pca.data$x), x = ~PC1, y = ~PC2,
  text = ~paste("ID: ", rownames(pca.data$x)),
  color = ~as.factor(deaths$dead),
  colors = c('#636EFA','#EF553B'), 
  type = 'scatter', 
  mode = 'markers') %>%
  layout(
    legend=list(title=list(text='color')),
    plot_bgcolor = "#e5ecf6",
    xaxis = list(
      title = sprintf("standardized PC1 (%.2f%% explained var.)", explained_variance_pct[1])),
    yaxis = list(
      title = sprintf("standardized PC2 (%.2f%% explained var.)", explained_variance_pct[2])))

for (i in seq(length(features))){
  fig <- fig %>%
    add_segments(x = 0, 
                 xend = loadings[i, 1], 
                 y = 0, 
                 yend = loadings[i, 2], 
                 line = list(color = 'black'),
                 inherit = FALSE, 
                 showlegend = FALSE) %>%
    add_annotations(x=loadings[i, 1], 
                    y=loadings[i, 2], 
                    ax = 0, 
                    ay = 0,
                    text = features[i], 
                    xanchor = 'center', 
                    yanchor= 'bottom')

}
fig

Task 10

Дайте содержательную интерпретацию PCA анализу. Почему использовать колонку ‘dead’ для выводов об ассоциации с выживаемостью некорректно?

To be done

Task 11

Приведите ваши данные к размерности в две колонки через UMAP. Сравните результаты отображения точек между алгоритмами PCA и UMAP

umap_comps <- data %>% 
  select(where(is.numeric), -c(ID, birth, year, exit)) %>%
  na.omit() %>% 
  recipe(~.) %>% 
  step_normalize(all_predictors()) %>%
  step_umap(all_predictors()) %>% 
  prep() %>% 
  juice()

combined_data <- cbind(umap_comps, pca.data$x[,1:2], deaths)
umap_plot <- combined_data %>%
  ggplot(aes(UMAP1, UMAP2)) +
  geom_point(aes(color = as.character(dead)), 
             alpha = 0.7, size = 2) +
  labs(color = NULL) +
  theme_bw()

pca_plot <- combined_data %>%
  ggplot(aes(PC1, PC2)) +
  geom_point(aes(color = as.character(dead)), 
             alpha = 0.7, size = 2) +
  labs(color = NULL) +
  theme_bw()

ggarrange(
  pca_plot, umap_plot, 
  nrow = 1, ncol = 2,
  widths = c(1, 1)
)

В целом, результаты получились похожими, хотя умершие субъекты лучше “кластеризуются” в UMAP.